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A CURVE FITTING METHOD FOR 
SOLVING THE FLUTTER EQUATION 



CHAPTER I 
INTRODUCTION 

In the optimization of a complex structure such as an airplane, 
one of the usual design constraints is that of flutter. Flutter is a 
hazardous phenomenon, for which the physics is not fully understood, 
in which an elastic body in an airstream, e.g. an airplane wing, 
becomes dynamically unstable. The first flutter investigation occurred 
about 1915, but an accurate theory of aircraft structure flutter was 
not formulated until 1935- Engineers in the recent past were able to 
design flutter free airplanes with. the existing established techniques, 
but these methods without the digital computer meant long cumbersome 
hand cal cul ati ons . As always time consumption meant money consumption, 
and with the age of new computers automated methods of flutter analysis 
were needed. Automated methods have now been applied to the different 
methods of flutter analysis, but these are basically iterative 
techniques that could be time consuming and therefore costly. 

It is the objective of the author, not to extend the theory of 
flutter, but to present a flutter velocity extraction technique which 
is relatively fast with respect to digital computational time and is 
very simple. The flutter analysis presented here is a part of a large 
optimization computer program that optimizes a simple three bay canti- 
levered box beam (Figure 9) for minimum mass subject to the flutter 
velocity constraint. The box beam was used since many airplane wings 
resemble a box beam structure. Two examples are presented. Case 1 is 



2 


the system resulting from what will be defined as the minimum variable 
parameters (P(l) - PMiN(l)). Case 2 is the system resulting from a 
set of variable parameters which are six times the minimum parameters 
(P ( I ) “ PMiN(l) x 6). These variable parameters (cross sectional areas, 
thicknesses, etc,) are listed in Table 1 1 1 of the Appendix. 

Before explaining this method of flutter analysis a better 
understanding of the overall problem of flutter should first be 
developed. The following chapter contains a discussion of flutter 
for a system with single degree of freedom followed by the broader 
concept of flutter for a system with many degrees of freedom. The 
formul at ion of. the basic f lutter equat ion is also presented in this 
chapter. The classical or established methods of flutter analysis, 
recent work in the area of flutter, and the statement of the problem 
are given in Chapter III. Chapter IV presents the derivatives of the 
eigenvalues ^ with respect to the reduced frequency k, which could 
have been derived directly from a paper by Rudisill and Bhatia (11), 
but which are given in detail since some of the complexity and under- 
standing might be taken from the problem. These derivatives make 
possible the computerized flutter technique explained in Chapter V. 

The results and conclusions of this technique are also given at the 
end of Chapter V. 
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CHAPTER I I 

FLUTTER AND THE FORMULATION OF THE FLUTTER EQUATION 

As defined in the Introduction, flutter is a hazardous phenomenon 
in which an elastic body In an airstream becomes dynamically unstable. 
The elastic bending of the wing is such that energy is added to the 
structure faster than it can be dissipated. The flutter speed and 
the frequency ^ are the lowest airspeed and the corresponding circular 
frequency of oscillation, respectively, at which a given structure 
flying under given atmospheric conditions will exhibit sustained simple 
harmonic oscillations (1, p. 5). Flight at velocities above the flutter 
airspeed will cause divergent oscillations or damage to the structure, 
while airspeeds below represent stable conditions. The flutter 
speed Is therefore a condition of neutral stability for the system. 

Theoretical flutter analysis many times consists of assuming in 
advance that all dependent variables are porportional to e ?a)t where <x> 
is a real frequency and i- (-] ) Combi nations of the velocity V and 
the frequency u for which this occurs must be found. This double com- 
bination thus leads to a complex or double eigenvalue problem, with 
two characteristic numbers determining the airspeed and the frequency. 
The simpler analogous situation is the free vibration of a linear 
structure in a vacuum, which is a real or single eigenvalue problem. 

In the past, flutter predictions have been made by many different 
methods. Some of the methods used were analog simulation, scaled 
dynamic models in wind tunnels and flight testing of full scale models 



by pilots. The last method requires that the airplane be flown at 
the predicted flutter velocity. The excitation of the structure at 
this speed sometimes causes divergent oscillations which could mean 
death to the pilot. 

Single Degree of Freedom 


The flutter problem discussed in this paper is one for a 
structure with many degrees of freedom; therefore, a simple system 
will first be explained. Bisplinghoff (2) considers the example of 
a rigid two dimensional airfoil in Figure 1. The airfoil is hinged 
at its leading edge so that it is elastically restrained from rotating 
about its leading edge due to the torsional spring with a spring con- 
stant equal to ft-lb/rad. The unstretched position of the spring 

corresponds to a zero angle of attack a. The equation of motion for 
this single degree of freedom system is 


where 


I a + K a 
a a 



2-1 


Ig-' - moment of inertia about the leading edge 
My - aerodynamic moment due to a(t) 
a = angle of attack and 
- spring constant. 

In order to produce flutter for this simple problem, some 
unrealistic assumptions must be made. These assumptions are that the 
air density p will be very small and the air foil will be heavily 
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weighted. First solve for the flutter condition by assuming as a 
solution for equation (2-1) the following, 

— i^t 

a = a^e 2-2 

where Oq is a constant angular amplitude displacement. Substituting 
equation (2-2) into equation (2-1) yields, 

I (a ^ i 2 aj 2 e ' a>t ) + K (a e la)t ) = M . 2-3 

a 0 a o y 

Upon dividing equation (2-3) by Trpb^o^a^e' ^ , where b is the reference 
semi chord, the result is 

* I w 2 

[j - {--£) ] + m = o 2-4 

u y 

. TTpD 

oc 2 

where (- — ) is the natural • frequency of torsional vibration in a 

a 1 i _> 

vacuum and m y = M y /irpb u a^e. represents the dimensionless aerodynamic 
coefficients. For a thin airfoil with small harmonic motion in two 
dimensional flow m y .may be written as 

m y = M a + l/2(L a + 1/2) + l/4(L h ) 2-5 

where L^, and are complex functions of the reduced frequency k. 
The quantity is therefore a complex function and can be split into 
real and imaginary parts, 


Re(m ) = 


tii 2 

U-r) - U 


lm(m ) ■ 0. ; 

y' < 

From equations (2-6) and (2-7) flutter occurs for the value of 
the reduced frequency k for which the out-of-phase component of the 



r" 

6 



K - spring constant of torsional spring 
• v t> = semi chord 

a - angle of attack 
V - velocity of airstream 


Figure 1. Rigid Two-Dimensional Airfoil 'Restrained about 
i ts Lead i ng Edge . 
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aerodynamic moment will vanish, provided that the 

corresponding in-phase part, Re(rn^), is of such magnitude that equation 
(2-6) y iel ds a pure real flutter frequency. The out-of-phase component 
is the only damping or instability for the airfoil. For Im(m^) 
less than zero or an m^ which lags the simple harmonic motion, energy 
is removed from the oscillations, resulting in damping of the airfoil. 
If Im(niy) is greater than zero then divergent oscillations result with 
energy added to the airfoil. 


A System of Many Degrees of Freedom 


By the finite element method a complex structure with many 
degrees of freedom is described as an assembly of structural elements 
connected by discrete or continuous attachments, and the displacements 
at nodes, including angular displacements, are represented as the 
actual displacements on the continuous structure. The interaction 
forces, including bending moments and torques between the discrete 
attachments, are represented as discrete joint forces. Figure 2 
shows how a cantilevered beam would be represented using the finite 
element method. The applied forces P., P_, the reaction force P , the 

I L f 

applied moments Mj, M^, and the reaction moment are replaced by 
discrete forces, H ]t Y ] , M 2> Y 2 , Y y M^, Y^, and noda) displacements 
V|f 0|, v 2 , e 2> v^, 0^, v^, 0^ at nodes (1), (2), and (3). 

For the system with many degrees of freedom and therefore many 
elements the elastic and inertia properties of these elements are 
determined in matrix form. It is from these properties that the 




Figure 2. Discrete Element Representa 
Beam. 
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matrix equations for static or dynamic equilibrium are derived. The 
motion of a discrete and linearly elastic system can be expressed in 
the form (l 1 ) 

W{U} + [K]{U>-{P} 2-8 

where, 

{U} = n x 1 vector representing nodal displacements for the 

assembled structure, referred to a global coordinate system, 
{U} = d {U}/dt = n x 1 vector representing nodal accel erat ? ons , 

[K] “ n x n stiffness matrix for the complete structure cor- 
responding to displacements {U}, 

{P} = n x 1 vector representing the total applied forces at the 
nodes,, corresponding to the nodes by {U}, and 
n « the number of degrees of freedom for the structure. 

As with the single degree of freedom system the theoretical flutter 
analysis for many degrees of freedom consists of finding combinations 
of the airspeed V and the circular frequency of oscillation cu for which 
the out-of-phase component of the imaginary part of the aerodynamic 
moment vanishes, producing simple harmonic motion. Based on the 
assumption of simple harmonic motion, the displacement vector {U} 
and the acceleration vector {U} may be expressed in terms of a constant 
amplitude vector {U^} as 

{U} = {U Q } e iut 2 _g 

and 

CU} = -w 2 {U Q } e iwt 
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The applied loads, P{t), may be expressed as 

P(t) — • oj 2 pb^ L [A]{U} 2-11 

where p, b, [A], and L are the air density, semichord, air force 
matrix, -and span length respectively. Substitution of equations 
(2-1G) and (2-11) into equation (2-8) yields, 

-w 2 e l<Jt [H]{U 0 } + [K] (U> = w 2 pb 4 L [A]{U> 2-12 

Equation (2-12) may be written as 

-w 2 [M]{U} + [K] {U} - oi 2 pb 4 L [A] {U} 2-13 

Finally equation (2-13) may be written as, 

[[K] -X([Hj + c[A]) ] (U) *-{0} 2-14 

4 2 

where c = pb L. and A = a> . Here X is referred to as the eigenvalue of 
equation (2-14). Equation (2--14) is the fundamental equation of 
f lutter. 

The sti.ffness matrix K and the inertia matrix M are functions of 
the variable parameters (cross sectional areas, thicknesses, diameters 

squared, etc.). The air force matrix A is a complex function of the 

\ 

reduced .frequency k, the air density, the mach number, and the semi- 
chord b. The reduced frequency k is defined by the relation, 

k - bw/V^ 

where is the flutter velocity. The matrices K and M are symmetric 
while matrix A is not. 

For a nontrivial solution of equation (2“l4) the flutter 
determinant | [K] - X( [H] +'c[A])| must be equal to zero. The known 
quantities of the determinant are matrices K and M which are determined 
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from the variable parameters of the structure and therefore the 
above determinant yields a polynomial in \ (=id ) with unknown complex 
coefficients which are in general functions of the reduced frequency 
k and mach number. The theory used for formulation of the air force 
matrix A for this paper is that of Theodorsen for subsonic incompress- 
ible flow. The use of this particular theory yields a polynomial from 
the flutter determinant which is a function of the reduced frequency 
k only. 

Real values of the reduced frequency k are assumed and the 
corresponding values of the eigenvalues are computed for which the 
flutter determinant is zero. A system with n degrees of freedom will 
have n distinct complex values of X corresponding to an assumed value 
of the reduced frequency k, if there are no repeated roots. Repeated 
roots are seldom encountered in a structural system. The flutter 
velocity analysis consists of finding values of the reduced frequency 
for which the imaginary part of the eigenvalues will equal zero. 'The 
flutter velocity is then computed from the following relation 

= bu>/k 2-15 

It is seen from equation (2-15) that although a value of 1/k 
may yield a real value for the eigenvalue X, it may not be the lowest 
flutter velocity, since is also a function of the circular frequency 
of oscillation w. The semichord b will be a constant. Thus, the 
lowest of the flutter velocities must be found for the system. In 
practice the lowest flutter velocity is usually found by checking only 



12 


the first three to five flutter velocities for Increasing values of 
1/k. No proof that this is always the case has ever been presented 
in the flutter literature, however. 



13 


CHAPTER III 
WORK BY OTHERS 

Lawrence and Jackson (5) discuss the flutter analysis methods 
classified as the American approach, the British approach, and the 
Richardson approach. In their paper they have shown the results of 
the three methods, compared them, and given their conclusions. It 
is the purpose of the author to discuss in this chapter what may be 
considered the classical techniques of flutter analysis and then 
briefly summarize the recent state of the art. The American and 
British techniques will be considered as the classical methods of 
analysts. The different methods of flutter analysis commonly used 
agree in regard to critical flutter airspeeds, provided the same basic 
data are used, but give different results for the stability decay 
rates at other airspeeds. 

British Approach 

In the United Kingdom the flutter equation is written in the 

form, 

M ^r- + t ? f D a ] + t D s ]] ^r + tV 2 [K a ] + [K S 1J- {UJ=0 3-1 

where 

[HI = st ructural and aerodynamic inertia matrix, 

{U} * eigenvectors or column matrix of generalized coordinates, 

V ~ a i rspeed 
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Vq =* reference airspeed, 

V = V/V Q , 

L » span length, 
t = V Q t/L, 

[D ] = aerodynamic damping matrix, 

[D ] “ structural damping matrix, 

[K ] - aerodynamic st i ffness matrix, 

Q . 

[K ] - structural stiffness matrix. 

A solution of the form {U} = {U} * 


constant amplitude vector and X , are the 

J 

of the form 


and 

is assumed where { U } is a 
resulting complex eigenvalues 


X . - y . + iuj. 
J J J 


where u is the exponential growth rate and ai is the circular frequency 
of oscillation. The system is stable if p is greater than zero and 
the state of flutter occurs when p = 0, or X = ito, or the eigenvalue 
is purely Imaginary. Equation ( 3*0 holds only for simple harmonic 
motion. 

- The aerodynamic stiffness matrix and the aerodynamic damping 

matrix D a are functions of the frequency v = w/V", and usually the mach 

number. The structural damping matrix is frequently set equal to 

zero since it is negligible compared with damping matrix D due to the 

a 

aerodynamic forces. Solutions of equation (3-1) are found by 

assuming a value of v for the calculation of matrices D and K and 

a a 

then solving for the ei genva 1 ues X from equation (3-1). The eigenvalues 
in general are not purely imaginary, and in general the computed value 
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of y - u/V, where o> comes from the imaginary part of the computed 
complex eigenvalue, do not line up or equal each other. When the 
computed eigenvalue is pure imaginary (p = 0) , then the assumed v 
and the computed v will line up. 

In order to find the lined-up values of v a graph is plotted of 
w obtained from the eigenvalues against. V, and the intersections of 
each curve wi th the line w = vV, where v is the assumed value of the 
frequency which was used originally to compute matrices D and K , 

3 d 

give the lined-up values of frequency and airspeed V. This method 
i s graph i ca 1 ly represented in Figure 3- The lowest critical velocity 
from the lined-up values will be the flutter velocity. 

The main disadvantage of this method is that the eigenvalues 
for a large range of assumed frequency v" values must be calculated 
before lined-up values can be found. Therefore, the British method 
will in general be time consuming and costly for application on the 
computer. 


American Approach 


The American method of analysis, commonly referred to as the 
V-g method, is quite different from that of the United Kingdom. The 
flutter equation is formulated as follows, 


i V [ D ] TT 2 . 

[[M] - —JL. - (I) [Kj - [K ]] 




OJ 


w 

where 1/v = V/w. It is seen from equation 
damping matrix is completely ignored and 


{U} = 

(3-2) 
9> a 


0 3-2 

that the structural 
fictitious structural 
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damping is introduced. The structural damping will actually be 
represented by ( g/ to) [ 3 and is supposedly just sufficient to maintain 
simple harmonic motion. The solution {U} = {U} e' WT is postulated 
leading to an eigenvalue problem to determine the complex eigenvalues. 

Once again an assumed value of V which is used for calculating 
the aerodynamic damping matrix D and the aerodynamic stiffness matrix 
yields in general the complex eigenvalues from equat ion (3"2). By 
separation of the real and imaginary parts of the eigenvalue X, g and 
oj may be obtained. The corresponding velocity V is obtained from the 
relation V = v/w. A negative value of g means that the system is 
damped and therefore, stable. The critical flutter airspeed occurs 
when g - 0. Figure 4 shows a graphical representation of the 
American method. It is seen from Figure 4 that the solutions for a 
particular value of oj/V 1 ie on a straight line through the origin of 
the u>-V plane. It is also seen that lined-up values of \7 play no 
part in this method as compared to the British method. Further com- 
parison with the British method shows that all computed values of the 
frequency v correspond exactly with the assumed values of V. Once 
again the American approach calls for an iterative method varying v* 
to find eigenvalues which will yield g = 0. From the American formu- 
lation of the flutter equation a plotting process to find zero values 
of g would be time consuming and costly, but the American approach 
does have the advantage of not having to find the lined-up values of 


v. 
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Recent Work 

Richardson's method is basically an extension of flutter theory 
to give correct subcritica? decay rates for a system with many 
degrees of freedom and multidimensional compressible aerodynamic 
forces. Richardson (8) formulates the general motion of a system 
as follows: 

2 T 

E [m3 - Em 1 3 3 + Ik s Hu} = v 2 [a(t - dt Q 3-3 

dt o o 

where the matrix A(t) is the indicial aerodynamic matrix, and is 

related to the aerodynamic damping matrix D , the aerodynamic stiff- 

a 

ness matrix K , and the aerodynamic inertia matrix M.. The indicial 
a 1 

aerodynamic matrix A(t) is approximated by a power series expansion 

which was a main point by Richardson in this approach. 

The solution of equation (3"3) by Richardson results in m 

coefficients of the matrices D and K , which represent the unsteady 

a a 

state of the growth of lift functions for an airplane travelling at 
speed V in air or the coefficients resulting from unsteady aerodynamic 
forces. The solution of equation (3~3) results in n(m + 1) roots 
instead of the usual n complex roots correspond ing to n degrees of 
freedom encountered in the American and British methods. The n 
roots represent the usual decayed oscillations and the nm roots are 
negative real roots representing exponential decays. The Richardson 
method, like the British and American methods, gives accurate critical 
airspeeds for the same basic data, but below the critical airspeed 
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the British and American methods are not as accurate as the Richardson 
method for computation of the decay rate. Of the three methods the 
American method must be considered the simplest to use for prediction 
of the flutter velocity. 

Others that have recently worked with flutter are Hassig (3) and 
Phoa (7)» Hassig applies a determinant iteration method to the British 
method of solution. In his formulation of the flutter equation he 
includes transfer functions for appl icat « on to hydraul i c controls and 
automatic controls systems. The method is an iterative automatic 
search procedure for the flutter velocity. Phoa uses Nyquisit's 
technique for his computerized technique. His approach shows the 
appi i cabi 1 i ty to systems which include feedback control systems. 
Although automated, the method is also an iterative procedure. 

Statement of the Problem 

Although automated techniques (Hassig, Phoa) have been applied 
to the various methods of flutter calculation, these methods are 
iterative procedures that could be time consuming and therefore, 
costly. Hassig's determinant method, however, must be considered as 
a more direct method of solution than the method used by Phoa. 

It Is the objective of the author, not to extend the theory of 
flutter, but to find a direct method of solution for the flutter 
velocity that is fast and simple, and apply this method to the 
optimization of a cantilevered box beam for minimum mass due to a 
flutter velocity constraint. 
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Rudisill and Bhatia (11) have derived the second derivatives of. 
the flutter velocity from the basic flutter equation (2-14). From 
their work and this formulation of the flutter equation (2-lA) the 
first and second derivatives of the eigenvalues with respect to the 
reduced frequency (dA^./dk, d Aj/dk ) can easily be derived. With 
these derivatives it Is the objective of the author to apply a curve 
fitting technique for which the flutter velocity could be found 
directly and simply. The expressions for the first and second 
derivatives of the eigenvalues with respect to the reduced frequency 
could simply be written down, but the author believes that some of 
the understanding and difficulty would be taken from the problem. 

s 

Therefore, t.hese expressions are derived in the following chapter. 
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CHAPTER IV 

FIRST AND SECOND DERIVATIVES OF THE EIGENVALUES 

In Chapter li it was seen that the fundamental flutter equation 
was written as, 

[[K] - X([Ml + c [A] ) ] { U } - 0 2-15 

For future derivation equation (2-15) will be written as, 

[[K] - A. ([M3 + [Aj)']{U.} = 0 4-1 

J J 1 

where the notation [A] has replaced c[A] since c ~ pb^L is a constant 
coefficient of matrix A. The subscript j indicates there will be n 
eigenvalues A^ and n corresponding eigenvectors U. for a single solu- 
tion of the flutter equation (4-1) where n is the order of matrices 
K, M, and A. 

From equation (4-1) It is seen that [[K] - A . ( [M] + [A] ) ] is 
actually a composite matrix of the stiffness, mass, and air force 
matrices. Transposing this matrix yields 

EEk] - Aj ( [m] + [a])] t - [[k] t -;.Aj([Mj T + Ca] t ) ] 4-2 

Since the stiffness and mass matrices are symmetric, 

[K] T • Ik] and [M] T = [M] . 

The resulting eigenvectors from the transposed matrix (4-2) will be 
denoted as V^.. The flutter equation may now be written as, 

EEk] - a j ( [m ] + Ea 3 t )]{v.} = o 

Transposing this equation gives, 

{Vj} T [[K] - A^ ([M] + [A])] = 0 2|-3 
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For a linearly elastic structure such as the box beam, the dis- 
placement at a node <5. can be related to the applied forces by the 
following relation, 


<$ . - C..F. + c. 0 F., + ... + c. F 
• 111 i 2 2 inn 


4-4 


The coefficients c.^ are known as the deflection influence coefficients 


In matrix notation equation (4-4) becomes 

{<5} = [C] { F} 


4-5 


Matrix C Is known as the matrix of deflection influence coefficients 
or is many times referred to as the overall structure flexibility 
matrix. If the inverse of matrix C exists then the applied forces may 
be found in terms of the displacements. In matrix equation form this 
becomes 


m = [c] _1 {6> j|-6 

If [C] exists then [C] ^ - [K] is the stiffness matrix. This 
relation is also true for the flutter equation (4-1). 

Premultiplying equation (4-1) by the flexibility matrix C yields, 
[[I] - yC]([H] + [A]) ] {lh } =0 4-7 

In order to get equation (4-7) in a form used for later differentiation 
define, 


Aj = 1/A. - 1/w 


4-8 


Multiplying equation (4-7) fay -X yields. 


t.-A [ I ] + [C] ( [M] + [A] ) ] {U. } = 0 
** J 


4-9 



25 


Rearranging the terms of equation (4-9) yields, 

C[C3([H] + [A]) - X. [I]]{U } = 0 4-10 

which will be the form of the flutter equation to be differentiated. 
For digital computational purposes this formulation of the flutter 
equation was found to be the fastest and the simplest. Thus, the 
foil owing equation develops from equation (4-10), 

{Vj } T [ [c] ( [m] + [A]) - T [I]} = 0 4-n 

With reference to equation (4-11) let 

[Fj] = [[C]([M] + [A]) - I [I}] 4-12 

From equation (4-12) the flutter equation may be written as 

[F ]{U } =0 4-13 

and 
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The flexibility matrix C and the inertia matrix M are functions 
of the variable parameters only and not the reduced frequency k. 

The air force matrix A and the eigenvalues X^ are however functions 
of the reduced frequency. Differentiation of equati on (4-12) with 
respect to* the reduced frequency k will therefore yield 

dtFjJ/dk = [C] [I] it- 18 


Substituting this result into equation (4-17), the followina 
result is obtained. 


{Vj} T [c] {u j} - { Vj> T -ii- [I] {U } = 0 it-19 

Solving equation (4-19) for the first de r ? va t i ve of the eigenvalues 

dX. 

with respect to the reduced frequency, — — L , gives 



{V.} T [C] {U.} 

_J dk J 

{V j } T [l]{U j > 


it- 20 


Since {V.} T [I]{U.} = {V.} T {U.}, 
d J J J 


equation (4-20) may be rewritten as 



{ V.} T [C] {U.} 

— J dk J 

{V J } T {U j } 


4-21 


The expression for the first derivative may further be simplified 
by making the associated row vector \l ^ orthonormal with respect to 
the eigenvectors U. or 


{V.} T {U.} = 1 • 
J J 


4-22 
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where (V.} T is orthonormal. The orthonorma 1 izat ion will be explained 
in the following manipulations. 

Consider the following equations, 

{Vj } T [[c] ( [M] + [A]) - X7 [I ] ] = 0 I»-23 

and 

[[C]([M] + [A]) - IT = 0 4-24 

where X^ and X. ate eigenvalues of the respective equations. Post - 
multiplying equation (4-23) by {U.} and premultiplying equation (4-24) 
by yields, 

{Vj} T [[C]([M] + [A]) - T [l]Hu.} = 0 • 4-25 

{Vj > T [ [C] ([M] + [A]) - X.[l]]{u.} = 0 4-26 

Subtracting equation (4-25 from equation ' (4-26) , 

{v j} T (Ij - VtiHV = o 4-27 

is obtained. From equation (4-27) i f X. and T. are distinct (i £ j) , 
then 

{V.} T {U.} =0 4-28 

J * ° 

then for the nontrivial solution 

{V. } T {U. } = d 4-29 

J J 

where d is a constant. Dividing both sides of equation (4-29) by d 
yields, 

{V.} T (U.} = 1 
J J 


4-30 
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where V, is the orthonormal associated row vector of { U « } with 
•* J 

respect to {IL}. Using this result and applying it to equation (4-21) 
the first derivative of the eigenvalues with respect to the reduced 
frequency becomes, 

dl. 


-L= fv/[C] ^{U > 


4-31 


Second Derivative of the Eigenvalue 
With Respect to the Reduced Frequency 

The starting point for the second derivative will be equation 
(4-15) which is, 

d[F.] d{U.} 

-dir { y + ty ik 1 - - ° <»-i5 

Differentiating equation (4-15) with respect to the reduced frequency 


d 2 [Fj 


dk 


d[F } d(U.} d [F . J d{U . } 

i— {u. } + J J— + _J J— . + r F } 

2 j dk dk dk dk H 


d 2 {U.} 


J— = 


dk 


= 0 4-32 


is obtained. Premul t iply ing equation (4-32) by’ {V.} T yields, 

_ T d[F.j d{U.) 

■* J L. 


- r d fF;.] _ T d[F.l d{U.} 

{V }' 1 {U.} + {V.} — 35-J -rp- L- + { v.} 

j d k 2 J j dk dk j' 


dk dk 


_ t rf2 { u -> 

+ {v riF ] A- - 0 

J J dk 

It was previously shown from equation (4-14) that 


4-33 


(v.r [f. 3 = o 

j j 


4-14 
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It is therefore also true that, 

<v/t F j] = 0 A-3A 

Grouping terms and using equation (4-34), equation (4-33) reduces to 
the fol lowing, 


{ V 


d 2 [Fj 


d[F.] d{U.} 


— f- {u.} + 2{v.r -nr L --^r ! - = o 

-1 j dk dk 


4-35 


Fj was defined as^ 


[Fj] = t[C] ([M] + £a]) - X [I]] 


4-12 


Differentiating Fj with respect to the reduced frequency k, yields. 


and 


d[F .] < r . -i ^A * 

i— = [c] - —J r 1 1 

dk LLJ dk dk i,J 


d [FJ j2r.i d 2 I. 

t~ = tC] ±iAi-„L [| 3 
dk dk Z 


4-36 


dk 


4-37 


Substituting d[Fj]/dk and d^[Fj]/dk^ into equation (4-35) yields, 


4-38 


Solving equation (4~38) for d A^/dk , the resulting expression is. 


d 2 X. 2[Vj> T [[C]^l - ^C]]^ + {V j} T [C]^M {U > 


dk' 


tVj} T [l]{U } 


4-39 


Since the denominator is unity the expression for the second derivative 
becomes , 
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d 2 A. 


dk 


— T HfAl ^ - d{U. } _ ,2 r I 

i= 2 { V } T [[c]^l- {7 > T [c]^M (u } 4-4o 

dk J 


The only unknown of equation (4-40) is d{lK}/dk, which may be 

written as a linear combination of vectors U, , 

h* 

d{U.} n 

— -7^— = E a., , {U, } 

& h=J jkh h 


4-41 


where a. kh is a scalar quantity and are eigenvectors of the n 
eigenvectors which form a complete set of vectors in a n vector 
space. Substitution of equation (4-4l) into equation (4-15) yields, 


d[f ] n 

^- { v + [F j ] h !, a jkh {u h j ■ ° ^ 

Premultiplying equation (4-42) by {V.} T where i ^ j yields, 

T d [ F . ] __ t n 

<V + ^ a Jkh (U h >- 0 4-4 3 

From equation (4-43) 

{V.} T [Fj] = {V.} T [[C]([Ml + [A])- I .[!]] 4-44 

is obtained. Postmultiplication by' {U^} yields, 

{V.} T [F J ]{U h } = {V.} T [[C]([M] + [A]) - I . [i3Hu h J 4-45 

From the basic flutter equation (4-10) the following relation may be 
written, 


[[C3([H] + [A] ) ] {U. } - X,[I]{U.} 4-46 

• J J 
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Rremu 1 t i pi i ca t i on by {Vj} of equation (4-46) yields, 


ru iT 


{Vj } ' [ [C] { [M] + [A] ) ] {Uj } = {V } T A [ I ]{U . } , 
but the quantity on the right side of equation (4-47) yields, 

{v.} T X.[i]{u.} = I. 

j j j j 

j 

from the orthonormality of {V^} with respect to (IK }. Thus, it is 
seen that, 


4-47 


rT7 -.T 


{VjJ'ttcKM + [a]) - a. [ may = (Vj } T ttc] ([m] + [a])]{u } 


- {v.fi.DKu f 


A-4'8 


or. 


{V.} T [F.]{U.} = a. - a. 

J « J J I 


4-49 


Applying the above relation to equation (4-43) gives 
- T dtF i :1 - _ 

{ V -*& {u j } + a jki (x i - y - 0 

y 

for h = i, since {V.} [Fj](U^} = C for h ^ i , Solving for yields, 


- { V. }T -^il {U j } 


d[F. ] 


a . i . 
jki 


X. - X. 
' J 


4-50 


where i £ j. The following relations are shown to give the case 

when i = j or the expression for a , even though these coefficients 

J K J 

will be seen to vanish in later derivations. 

It was shown previously in equation (4-30) that {7.} T (U.} = 1. 

It can also be shown in a similar manner that { Uj } can be normalized 


so that 
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and that 


{U.} T {U.} 1 , 

J J 


{U.} {U, } = constant * b., . 
J h jh 


4-51 


4-52 


Differentiating equation (4-50 with respect to the reduced frequency 
k y iel ds , 

lT T dtU.J 

4-53 


_ d{U.} 

5k {u j } + { y — ah- 0 


From equation (4-40 it was seen that d{Uj}/dk could be written as 
linear combination of vectors {U^}, thus the following relation for 
d{U.} T /dk is, 


d{U.} n 

= a jw <V T 


4-54 


Substituting equations (4-4l) and (4-54) into equation (4-53) yields, 

4-55 

Factoring out the a ^ coefficients from this equation the result is 

n T 

U.) + 2 a.. . (11 0(11 > + 7 , 

jkj 

Since {U^} {UO = b^. and {U^}"^{U^} = b^, the expression for a.^. is, 


T n T 

* a ikh {U h } {U ! } + 2 a... = 0 

h=l J Kn h J h= ] Jkh j h 


h f, a jkh {U h }T{U J } + j, -j kh «Wj> T f*h> + 2 ^: = 0 '•-56 

h^j h^j 


J jkj ~ ~ 2^ } a jkh^ b hj + b jh^ 

h^j 


4-57 
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The second derivative of the eigenvalues with respect to the 
reduced frequency was expressed as, 


a 2 r. 


dk 


* t , r a i dX , d {U . } , 2 r 

i-= «V {v j } tc i^r- { Uj } ^ 0 

d{U,} 


Let us focus attention to the terms containing — tt— of equation 

('.-'to). 


T , r ., dX. d{U.) T d(F.] d{U.) 

2 { Vj} T t[c]^I - ^ [ 1 ]]^-= 2 tv.} T 


4-62 


Substitution of the expression for the f i rst derivative of the 
eigenvalues with respect to the reduced frequency into equation (4-62) 
yields, 

2 <V -dk^ sb = 2 {V.} T UC]^i - (V.> T [C]^I (U.H,]]-^- 4-63 

where {V.}^[C] {U.} is a scalar. 

J dk j 

Since {Vj}^[c3- ^ {IL} is a scalar, equation (4-63) may be written as 
— - T^^i^ i ^ — j d [ Al d{U i* — T 

2{v j } V'-dk 1 -- 2{V j } [c3 ^T 1 “ W 


j 

it was shown previously that {V^} {Uj } = 1 from equation (4-30). 
Taking the transpose of equation (4-30) yields 

{U.} T {V.}=1 
J J 

It is also, therefore, true that, 


J J 


dk 


4-65 


{V.} T — J- = -- J' {V.} 
j dk dk j 


d{U.} d{U .} 


4-66 
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d{U.} 

With this result ^ may be written as 
d{U.} n 

- 2 a , CU,} + a... {U.} 


dk 


h :-jkh jkj 

h^j 


4-58 


or substituting the expressions for ^ and a^j yields, 

d [F . } 


dk 


d{U.} n {V, } T — {u.HU. } , n 

■ J -~- E — h - ' dk _ - 1 h - Us a... (b. . + fa )){U } A- 59 

h=l .. X. - X. 2 h=l Jkh hj Jh J 

hi 4 j J hsfj 

d{U,} T 


In a similar manner — — — may be written as, 


- T d[F i ] T T 

d{u.> T ,fVV‘V { V 


dk 


J_= - 


h=l 

hj*j 


-l 

A, ~ X. 
h J 


- k s a .. .(b. . + b ..) ) (U. } T . 4-60 

2 h=l Jkh hj Jh J 
h/j 


Postmul t iply ing equation (4-59) by (Uj > T and premu] t iply ing equation 
(4-60) by {LL} and subtracting the two equations the result is, 


T - T d[F i ] T 

^ {u } T ■ [ ;.- ] dtu i> . 2 { y r V-<V<W T 

dk tu j } { V dk h 2 , x - 

wj h j 


4-61 


+ r 
h=1 
h^j 


- T d[F | ] T T 

n (V -sr 1 - {o.}{u.} T {o. r 

h dk j i h 


A. - A. 
h j 


From equation (4-61) it is seen that the coefficients a.. , have 

Jkj 


vanished from further manipulat ions . 
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Substitution of (U } T {V.} for "1 M and the expression for {V.} T d{U,}/dk, 

•J J J J 

—Equation (4-64) becomes, — 

_ _d[F ] d(U } __ f . d{U. } _ d{U. } T 

2{v j ) ^r L -V- = 2{v j } ^-i 1 hnr^V <V - ^ J >-dir i — {? j i} 

4-67 


Factoring out {Vj}, equation (4-67) may be written as, 

__ T d[F.] d{U . } . d{U . } _ d{U.} T 

2{ y 2{v j } [c ^~ {- dir L ' { v ■ {u , } ~dif- Hv i } 


68 


r " ^ d{U, } d{U. } T 

Equation (4-61) gives the expression for { ^ — {U } - {y } — — J } 

dk j j dk 

Substitution of this equation into equation (4-68) yields, 

217 j }T ~dI^ V- - 2{V j {- E — h . } J HU h HU j }T 

h=) A, - A. 

h^j h j 

r^F.] T T 

t : { v - 3 ^ < - 

+ E U J Q_ }{v } 

h=l A- A. J 

h^j h J 


4-69 


It was previously shown that {V.} T {U.} = 0, where i * j. Thus 


{iy'{V } = 0 


4-70 


Where h ^ j. . Equation (4*70) causes the second term of equation (4-69) 
to vanish. The resulting equation is, 


h-1 

h^j 


A, - A. 

h J 


4-71 
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From equation (4-18) it can be shown that, 

_ T d[F ] ... _ . m __ 

{ V wftv 

From the result of equation (4-72) , equat ion (4-71 ) may be written 
as. 


4-72 


:[?J T d[F i ] „ 2 S ({V i >T[C]j # f »h»«V h } T [C]^I-{U,)) 


j dk dk 


h- 1 

h*j 


A. - A, 
J h 


4-73 


Substi tutlon of equation (4-73) into equation (4-40) yields thi 
expression for the second derivative of the eigenvalue with respect 
to the reduced frequency as, 


d 2 A. 

J 

<,k 2 


. ,7.,'tciilil »» J > T lc]|M (Uh , 1 „ ?| ,, T [ cidM,„,,| 


d[A] 


dk 


h-1 


X. - A, 
J h 


4-74 


The following chapter explains the author's method of finding the 
'Critical flutter velocity. Without the first and second derivatives 
of the eigenvalues with respect to the reduced frequency the method 
used would be impossible. - 
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CHAPTER V 

A COMPUTERIZED CURVE FITTING FLUTTER ANALYSIS 


From Chapter IV it was seen that the flutter equation could be 
written as, 


Uc]([M] + [A]) - I.UlKo.} = 0 


If- 10 


where C, M, and A were defined as the flexibility, inertia, and air 
force matrices, respectively. {Uj} were defined as the eigenvectors 
or a column vector of generalized coordinates. The solution of 
equation (4-10) for an assumed value of the reduced frequency k 
yielded j complex eigenvalues of the form 


A. = 1/A. « x- + i 
J J J J 


5-1 


1 . 1/2 

where Xj “ * > - ("D » and represents the imaginary part of 

. “j J 

the eigenvalue. The values of are labeled in descending order of 
Xj. Each j subscript corresponds to a mode of oscillation. 

The eigenvalues Xj are a function of the reduced frequency k, 
which is a function of the semichord b, the circular frequency of 
oscillation oi, and the velocity V by the relation. 


The eigenvalues can then be said to be a function of b/k or V/u). If 
the imaginary part of the eigenvalues #j are p.otted against v. , 
where Vj = V ./ajj a set of typical curves might like the curves of 
Figures 5 and 6. Figure 7 gives a representation of a set of typical 

Xj curves. It is seen from Figures 5 and 6 that modes 2, 4, and 6 
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produce ijn curves that cross the v axis. At these points where A_ 
is pure real or is zero, flutter will occur. The lowest velocity 
corresponding to a pure real X^ will be the critical flutter velocity 
for the system, where the flutter velocity Is computed from the 
relation, 


V f = (Xjf'V 


5-2 


From a structural synthesis point of view it is computationally 

inefficient to make plots of versus v. for each flutter analysis. 

J J 

Instead of using the plotting procedure, the crossings of ^ may be 

approximated by fitting cubic or quadratic equations to iJj. from com- 

2 2 

puted values of if/ j , di|/j/dv, and d tyj/d v. The first and second 
derivatives of the eigenvalues with respect to the reduced frequency 
were derived in Chapter IV. in order to use the curve fitting 
technique, the first and second derivatives of the eigenvalues with 
respect to v must be derived. 

Let D = b/k, where b is the reference semichord and k is the 
reduced frequency. Then, 

... . . dD = dk 5-3 

k 

or. 



5-4 





i 



Figure 6. Imaginary Part of the Eigenvalues, Modes 4, 5, 6. 


o 




L og( x .xlO 
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The first derivative of the eigenvalue with respect to b/k may be 
expressed as, 

dX. dX. 




dk 


dk , b 
d k 


5-5 


where. 


S i = {Vj} T [C]^fI {Uj} 


4-31 


Since b/k = V/oj = v, the expression for the f i rst derivative of the 
eigenvalue with respect to \> may be written as, 

^ <y . « 

Differentiating equation (5-5) with respect to b/k yields, 

.^ii_ - ill a, 2 + 5.a_ 

d(^) 2 . dv 2 dk 2 d(£) dk d(£) 2 

where dk/d(b/k) = -k /b and d k/d(b/k) — 2k^/b^. Substitution of the 

relations for dX./dk and d^A./dk^ into equation (5-7.) yields, 

J J 

... — £ = (V. } T [e]^Ii3- {u. } ■ 

dv Z b 2 J dk J 


+ (-£-) 2 {V } T [C]^M {U.} 
b J dk 2 J 

rTT \Tt>i_d[A] 


+ 2Z " { V i }T [Cl^-{U h MV h > T tCl^l { U |} 

h- ) A. - A, 

h^j 


5-8 
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Curve Fitting Technique 


-The~ipj curves are approximated by cubic equations of the form, 

= A.v + B.v^ + C.v^ 5*9 

J J J J 

The graph of cubic equations (5“9) pass through the origin and the 

coefficients A., B., and C. are evaluated such that the curves of 
J J J 

equation (5*9) pass through the points ( v 0 »^j) and (v, tj/.). The first 
derivative of with respect to v is, 


d^. 

* A. + 2B.v + 3C.v' 
dv j j j 


5*10 


The curves of equation (5*10) pass through the point (v q , d^?/dv) 

where i^? and d\p?/dv denote values of and d^/dv evaluated at 

v = v respectively, and ip. denotes evaluated at v = v /M, where 
o J J o 

M Is some large number, say 10 or 100. is an assumed trial value 
of v. 

The coefficients A., B., and C. are calculated from the 

J J J 

following equations, 


= A.v + B.v ^ + C.v ^ 

J JO jo jo 


d$?/dv =A. +2B.v + 3C .v ^ 

J J JO JO 


5*11 

5*12 


and, 


— — _2 —3 

ifc. = A.v + B .v + C ,v 
J J J J 


5-13 



Using Cramer's rule the coefficients of equation (5~9) may be 
expressed in the forms, 

-2 2x . . o , . /— 2 3 - 3' 


r , .v 4 + <|/?(2v 3 - 
I o r j 

,r ..3 ^ .,0,,- 2 —3- 


Aj - + ^j(2v^v^ - 3v^v^) + d^?/dv(v^v^ - vv^) }/D 


Bj - {-24»jV^ + (3vv^ - v 3 ) + dTj/?/dv(v^v^ - vv^) }/D 

C. = + *j>.(v - 2v\) ) + di//?/dv(vv ^ - v^v ) }/D 

J J ° J o T J o o 


5-14 

5-15 

5-16 


and, 


- — 4 —2 3-32 

D - vv - 2v v + v J \) 
o o 


5-17 


Since is zero at v - = 0, then the other two crossings of ipj may be 

-found by setting ^ in equation (5"9) equal to zero, thus 


A. +B.V+C.V =0 
J J J 


5-18 


From the quadratic formula the roots may be expressed as, 

v = {-Bj ± (Bj - i(AjCj) ,/ 2 .}/ 2 Cj 


5-19 


2 

provided is not zero and (B^. - bA.C^) is greater than or equal to 
zero. If v from equation (5~ 1 9 ) is not pure real and positive then 
curve fitting by the use of a quadratic equation of the same form as 
equation (5“l8) may prove successful. 

\ 

The coeff icients (Aj, B ^ , Cj) of the quadratic equation, 


“ A. + B.v + C.v‘ 
J J J J 


5-20 


are computed for curves which pass through points (v , it/?) and 

° J d*? 

have first and second derivatives which pass through points (v — -J) 

o*dv 
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2 2 

and (v q , d i^/dv ), respectively, may be found from the following 


equations , 


= A. + B ,v + C .v ^ 
J J jo jo 


dijj./dv = B. + 2C,v and, 

J J Jo 


5-21 

5-22 


2 o 2 
d i|>./dv = 2C . 

J J 


5-23 

From equations (5-21) through ( 5“ 23) the coefficients may be expressed 

as . 


and, 


2C . = d 2 tp?/dv 2 
J J 


B . - ./dv - 2C ,v • 
J J jo 


A . ~ ~ B . v - C . v ^ 

J J JO JO 


5-2V 

5-25 


5-26 


Substitution of these coefficients into equation (5-19) yields the 
roots of equation (5"20). 

If the roots of a mode are not real and positive, then the curve 

.fit for that mode fails. If one of the C. for the cubic equation is 

J 

zero then the approximate crossing for that mode may be found by using 
the Newton-Rapbson 1 s method for finding the root, provided d^?/dv is 
not also zero, then, 


v = v - ^?(d\p?/dv) 
0 J J 


5-27 


It was noted previously that the flutter velocity is computed 
from the relation. 


- 1/2 

s } y j 


a (x ; ) ' v. 


5-2 



Since the flutter velocity is computed from the predicted root v j , 

then the flutter velocity will be correct only if the corresponding 

2 

Xj value is used. The values of Xj = lAu at the crossings may be 
estimated from the first three"terms~crf" the Taylor expansion of x j > 
then 


Xj - Xj + (v-v o )dx?/dv + 0.5(v-v o ) 2 d 2 Xj/dv 2 5-28 

From the cubic and quadratic fit only the positive roots which also 
have positive slope di^./dv at the crossings, will approximate the 
lowest flutter velocity. 

A simplified flow diagram for the curve fitting flutter analysis 
is shown in Figure 8. Following Figure 8 it is seen that an initial 
value of v = v q , M, and some tolerance e are assigned. The quantity 
v - v^/M is then calculated which will cause the cubic fit to have a 
negative slope near the origin. The imaginary parts of the eigenvalues 
are then computed at “ for all n modes. = 10^° and v* = 0 are 
the upper and lower limits of the flutter velocity and the v value, 
respectively. it should be noted that each mode of oscillation n is 
checked for a crossing of i|k and a critical velocity is computed for 
each crossing. The velocities which are computed from the predicted 
crossings and the approximate value of x* (from Taylor's expansion) 
are compared to find the smallest velocity, which is then set as the 
new upper limit V f . The flutter velocity V f predicted from the 
initial assumed value is then compared with the velocity 
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If the absolute value of (V^ - V )/V is greater than some small 

t o o 

-number e, then the computations are repeated with set equal to 
the predicted crossover value of v which corresponds to the 
lowest critical velocity computed. If a positive crossing with 

it 

positive slope di/^/dv is not predicted then v ~ 0 wi 1 1 terminate 
computations. Also if the number of flutter iterations, for which 
there is no convergence, exceeds some number (e.g. 5), then computa- 
tions are terminated. 

Examples 

The data for equation (4-10) were generated for the cantilevered 
box beam (Figure 9) where the' design parameters were those listed in 
Table III. A digital computer program was written which implemented 
the flutter velocity solution previously explained. 

Case 1 will be considered as the minimum set of parameters (P(i) = 
PMlN(l)) and Case 2 as the minimum set of parameters times six 
(P(l) - PMIN(l) x 6). For the minimum set of parameters an initial 
value of v q = 1 was assumed and the program found a flutter velocity 
of 715.6 feet per second at v - 7-907 feet per radian, and for e 
equal to 0.05. The computations of the flutter velocity were repeated 
for values of ~ 2, 3, 4, ..., 20. In each case the program con- 
verged to a flutter velocity between 715.6 and 716.3 feet per second. 
The total time for the twenty flutter velocity calculations was 
245.53 seconds or an average of 12.276 seconds per flutter velocity 
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computation. From Table I it is seen that the maximum number of 
iterations to find the flutter velocity is five... For initial values 
of V o ec l ua ^ 7.0 and 8.0 the algorithm converged to an answer in 
two and one iterations respectively. 

The flutter velocity analyses were repeated for design parameters 
of Case 2 In a similar manner. The total time for twenty flutter 
calculations was 210.83 seconds, or an average of 10.5** seconds per 
-flutter velocity analysis. From Table II it is seen that for 
assumed values of 9, 10, and II that the algorithm converged to 
■the answer on the first attempt. The average time for a flutter 
-velocity analysis iteration for Casesl and 2 was 3.91 seconds. All 
computations were performed by an IBM 360, Model 50 digital computer.. 

Figures 12-35 show the actual computed ij/j versus v curves and 
the approximation of each curve as generated from the algorithm. 

The assumed values of v q - 5.0 and = 9.0 were used for Case 1 and 
Case 2 respectively. These values (v q ) were chosen since they lay 
between the crossings of different modes. Also, individual graphs 
are shown since in many cases the cubic or quadratic fit would lie on 
top of the computed curve for the majority of v values given. 

Tables V-VII give the numerical values of the plotted curves and 
the computed airspeed V for the reader to see a numerical comparison 
between the computed curve and the approximation. . 
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Table 

1. Convergence for Twenty Assumed Values 

v for Case L. 

0 

of 

Assumed 
Va 1 ue 
V 

o 

Predicted 
ROQt 
V ' 

Pred i cted 
Veloci ty 

V f 

Predicted 

Mode 

n 

1.0 

9.278 

869.2 

2 


7.762 

701.5 

3 


... 7.907 

715.6 

2 

2.0 

6.097 

554.2 

2 


7.416 

659.0 

2 


7.91 1 

714.8 

2 


7.906 

715.6 

2 

3.0 

5.641 

500.7 

2 


7.089 

627.0 

2 


7.885 

710.2 

2 


7.906 

715.6 

2 

4.0 

5.952 

523.9 

2 


7.326 

649.7 

2 


7.911 

714.4 

2 


7.906 

715.6 

2 

5.0 

6.611 

582.3 

2 


7.742 

692.5 

2 


7.908 

715.7 

2 

6.0 

7.374 

654.1 

2 


7.916 715.1 2 

7.906 715.6 


2 
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Assumed 

Value 

V 

o 

Predicted 

Root 

v" 

P red i cted 
Vel oci ty 

v f 

Predicted 

Mode 

n 

7.0 

7.886 

709.3 

2 


7.906 

715.6 

2 

8.0 

7.907 

715.6 

2 

9.0 

9.000 

1328.0 

3 


8.006 

724.9 

2 


7.907 

715.6 

2 

10. 0 

7.596 

689.0 

3 


7.915 

716.0 

2 

11.0 

7.268 

664.5 

3 


7.928 

715.5 

2 


7.906 

715.6 

2 


7.080 

653.7 

3 


7.926 

713.7 

2 


7.906 

715.6 

2 

13.0 

7.060 

658.7 

3 


7.931 

723.9 

2 


7.906 

715.6 

2 

14.0 

7.168 

675.4 

3 


7.941 

715.9 

2 


7.906 

715.6 

2 

15.0 

7.367 

700.3 

3 


7.937 

717.0 

2 


7.906 

715.6 

2 
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Assumed 

Value 

V 

o 

Predicted 

Root 

v* 

Pred i cted 
Ve loc i ty 

V f 

16.0 

7.628 

730.9 


7.917 

716.3 

17.0 

4.827 

471.7 


6. 569 

575.9 


7.861 

701.5 


7.906 

715.6 

18.0 

5.009 

494.3 


*6.725 

589.6 


7.933 

' 709.8 


7.906 

715.6 

19.0 

5.209 

518.8 


6.094 

605.4 


7.979 

716.2 


7.907 

715.6 

20.0 

5.422 

544.8 


7.099 

623.1 


7.987 

719.2 


7.907 

715.7 


. Predicted 
Mode 
n 

3 
2 

4 
2 
2 
2 
4 
2 
2 
2 
4 
2 
2 
2 
4 
2 
2 


2 
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Table If. Convergence 

v for Case 
0 

for Twenty Assumed 
2. 

Values of 

Assumed 

Pred i cted 

Predicted 

Pred ? cted 

Value 

Rogt 

Veloc i ty 

Mode 

V 

o 

v ' 

V f 

n 

i.O 

33.17 

9836.0 

5 


4.39 

352.4 

2 


13.63 

1180.0 

2 


8.99 

793.7 

2 


10.15 

870.4 

2 

2.0 

7.28 

2633.0 

5 


10.34 

892.2 

2 


10.16 

870.1 

2 

3.0 

28.96 

2126.0 

2 


4.28 

373.3 

2 


13.74 

H95.0 

2 


8.91 

787.7 

2 


10.15 

870.5 

2 

4.0 

14.63 

1270.0 

2 


8.04 

755.2 

2 


10. 19 

876.7 

2 


10.16 

870.1 

2 

5.0 

II .94 

1049.0 

2 


9.95 

857.3 

2 


870.1 


10.16 


2 
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Assumed 
Va 1 ue 

V 

o 

Pred Icted 
Root 

V* 

Predicted 
Velocl ty 

v f 

Pred Icted 
Mode 
n 

6.0 

10.87 

948.9 

2 


10.14 

868.8 

2 

7.0 

10.38 

899.2 

2 


10.16 

870.0 

2 

8.0 

10.18 

876.9 

2 


10.16 

870.1 

2 

9.0 

10.14 

870.2 

2 

10.0 

10.16 

870.1 

2 

11.0 

10.12 

868.2 

2 

12.0 

9.93 

856.0 

2 


10.16 - 

870 .O 

2 

13.0 

9.47 

826.0 

2 


10.15 

869.8 

2 

14.0 

8.75 

774.7 

2 


10.14 

870.8 

2 


10.16 

870.1 

2 

15.0 

7.89 

708.9 

2 


10.17 

877.7 

2 


10.16 

870.1 

2 

16.0 

7.12 

644.0 

2 


10.28 892.9 

870 . ] 


10.16 


2 



Assumed 

Value 

V 

o 

Predicted 

Root 

V'* 

Predicted 
Veloci ty 

v f 

P red i c ted 
Mode 
n 

17.0 

6.57 

591.9 

2 


10 . bl 

911.5 

2 


10.16 

870.0 

2 

18.0 

6.21 

555.2 

2 


10 . 

927.6 

2 


10.15 

869.8 

2 

19.0 

6.00 

531.4 

2 


10.62 

938.3 

2 


10.15 

869 .6 

2 

20.0 

5.90 

517.4 

2 


10.65 • 

943.6 

2 


10.15 

869.5 

2 
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Cone? us i ons 

For the two examples presented here the flutter velocity 

analysts program always converged to an answer for values of v chosen 

o 

between 1 and 20 and did not show any tendency to diverge from the 
solution, however, this does not guarantee that the method will be 
reliable for all systems. It may be possible that some other type of 
fit besides a cubfc or quadratic equation would give more accurate 
results, but the former was used because of simplicity. The value of 
M = 10 was used for the cases presented, but for some systems this 
number may be, either too large or too small and might yield a negative 
slope from the origin. 

It should be noted that the second derivatives of ^ are not 
needed unless the cubic equation fails to give a real positive value 
of v, however, the computer program computed all of the second 
derivatives even if they were not needed. Some computati onal time 
may be saved by computing only those values of d 2 ^/dv 2 which are 
needed. For the cases shewn the cubic was used for four of the 
six modes. 

One problem that did not effect convergence for the cases 
explained was that of “switching of the modes" for the ip. values. 

The subroutine which computes the eigenvalues arranges the eigenvalues 
Xj in descending magnitude of the real part of the eigenvalue x . . 

This is done since the Xj values for a particular mode vary only 
slightly as can be seen from Figure 10 or 11. The imaginary part of 
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the eigenvalues ^ . , however, for a mode may change from point to 

point by a considerable amount. The ordering of the eigenvalues by 

X-> therefore, may cause the ip’ value computed from v not to 
J J o 

correspond with the ^ value computed from v = v thus resulting in a 
bad fit for the actual shape of that mode. From Figure 11 switching 
is seen to occur for values of v = 3.0 and v = 14.0. 

it is the belief of the author that the flutter technique 
presented here is a basically simple and efficient method for finding 
the flutter velocity. For the designer with a feeling for the range 
of values to try, this method should be quite successful. 
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APPENDIX A 


Box Beam and Variable Parameters 
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Table III. Variable Parameters of Box Beam. 


Case 1, P ( I ) = PMIN(I), 


Bay 

No. 

Areas of 
Longi tud ina 1 s 

(ft 2 ) 

Front and Back 
Web 

Thi cknesses 
(ft) 

Top and 
Bottom Skin 
Thicknesses 
(ft) 

Rib 

Thi cknesses 
(ft) . 

1 

0.002315 

0.001111 

0.0005555 

0.0005555 

2 

0.002315 

0.001111 

0.0005555 

0.0005555 

3 

0.002315 

0.001111 

0.0005555 

0.0005555 

Case 

2 , P(!) = PMlN(l) 

x 6. 



Bay 

No. 

Areas of 
Long! tudi na Is 

(ft 2 ) 

Front and Back 
Web 

Th i cknesses 
(ft) 

Top and 
Bottom Skin 
Th i cknesses 
(ft) 

Rib 

Thicknesses 

(ft) 

1 

0.01389 

0.006667 

0.003333 

0.003333 

2 

0.01389 

0.006667 

0.003333 

0.003333 

3 

0.01389 

0.006667 

0.003333 

0.003333 



APPENDIX B 


Plots and Corresponding Tables of Computed 
Curves and Cubic or Quadratic Fit 



2.yO 3.20 





Table IV. Real Part of the Eigenvalues 
for Case 1. 


V 

ft/rad 

MODE 1 
(xl03) 

MODE, 2 
(xIO 4 ) 

1.0 

M3 

0.50, 

2.0 

1.09 

0.62 

3.0 

1.05 

0.78 

4.0 

1.01 

0.96 

5.0 

0.98 

Ml 

6.0 

0.97 

M9 

7.0 

0.98 

I.23 

8.0 

1.01 

1.22 

9.0 

1.05 

0.99 

10.0 

1.11 

M7 

11.0 

M9 

1.15 

12.0 

1.29 

1.14 

13.0 

1.49 

1.14 

14.0 

1.54 

1.14 

15.0 

1.69 

1.14 

16.0 

1.86 

1.14 


MODE 3 
(xlo5) 

0.38 

0.37 

0.37 

0.40 

0,46 

0.55 

0.68 

0,86 

1.19 

1.32 

1.57 

I.85 

2.15 

2.47 

2.82 


3.19 
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Table IV. Real Part of the Eigenvalues 
for Case 1. 


V 

ft/ rad 

MODE 4 
<x10 5 ) 

MODE 5 
(xl05) 

MODE 6 
(xlO 6 ) 

1.0 

0.76 

0.62 

2.96 

2.0 

0.86 

0.70 

3.35 

3.0 

0.97 

0.85 

3.50 

4.0 

1.17 

0.98 

3.44 

5.0 

1.47 

1.04 

3.31 

6.0 

1.86 

1.07 

3.19 

7.0 , 

2.35 

1 .08 

3.07 

8.0 

2.92 

1.09 

2.96 

9.0 

3.58 

1.09 

2.86 

10.0 

4.34 

1.10. 

2.79 

11.0 

5.19 

1.10 

2.72 

12.0 

6.14 

1.10 

2.66 

13.0 

7.18 

1.11 

2.61 

14.0 

8.32 

1.11 

2.57 

15.0 

9.57 

1.11 

2.53 


16.0 


10.90 


1.12 


2.50 
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Table V. Real Part of the Eigenvalues 
for Case 2. 


V 

ft/rad 


MODE. 1 
(xlO^) 

MODE. 2 
(xIO^) 

MODE 3 
(xl05) 

, 1-0 


8.91 

0.39 

3.07 

2.0 


8.85 

0.41 

3.04 

3.0 


8.76 

0.46 

3.01 

4.0 

_ 

8.66 

0.52 

2.98 

5.0 

- 

8.55 

0.61 

2.94 

6.0 


8.45 

0.72 

2.92 

7.0 ^ 


8.34 

0.85 

2.90 

8.0 


8.25 

1.00 

2.93 

9.0 


8.17 

1.17 

3.02 

10.0 


8.11 

1.34 

3-17 

11.0 


8.09 

1.50 

3-37 

12.0 


8.10 

1.66 

3.63 

13.0 


8. 16 

1.81 

3.92 

14.0 


8.26 

1.94 

4.26 

15.0 


8.41 

2.06 

4.64 

16.0 


8.59 

2.18 

5.06 



Table V. 


Real Part of the Eigenvalues 
for Case 2. 


V 

ft/rad 

MODE 4 
(xlO 5 ) 

MODE 5 
(xlO 5 ) 

MODE, 6 
(xiO 6 ) 

1.0 

0.49 

0.58 

2.36 

2.0 

0.53 

0.58 

2.51 

3.0 

0.59 

0.57 

2.77 

4.0 

0.67 

0.56 

3.16 

5.0 

0.77 

0.57 

3.59, 

6.0 

0.91 

0.61 

3.88 

7.0 

1.06 - 

0.67 

■4.06 

8.0 

1.20 

0.76 

4.19 

9.0 

1.32 

0.86 

4.29 

10.0 

1.41 

0.98 

4.39 

n.o 

1.49 

1.12 

4.47 

12.0 

1.56 

1.28 

4.54 

13.0 

1.61 

1.46 

4.60 

14.0 

1 .66 

1.67 

4.65 

15.0 

1.72 

1.87 

4.69 

16.0 

1.77 

2.10 

4.72 



MODE 1 (XIO 














Table VI. Computed i{j. and Quadratic Fit for 


Case 1 and s> 


' V 
ft/rad 

COMPUTED 

CURVE 

<xl0 3 ) 

0.0 

0.0 

1.0 

-0.15 

2.0 

-0.32 

3.0 

-0.51 

4.0 

-0.73 

5.0 

-0.97 

6.0 

-1.23 

7.0 

-1.50 

8.0 

-1 .73 

9.0 

-2.07 

10.0 

-2.36 

n.o 

-2.66 

12.0 

-2.96 

13.0 

-3.27 

14.0 

-3.58 

15.0 

-3.90 

16.0 

-4.22 


QUADRATIC 

v o = 5 *° 
(xlO 3 ) 

COMPUTED 

AIRSPEED 

ft/sec 

+0.09 

0.0 

-0.09 

29.7 

-0.29 

60.6 

-0.50 

92.6 

-0,73 

126.0 

-0.97 

160.0 

-1.23 

192.0 

-1.50 

223.0 

-1.79 

252.0 

-2.10 

278.0 

-2.42 

300.0 

-2.75 

319.0 

-3.10 

334.0 

-3.47 

346.0 

-3.85 

357.0 

-4.25 

365.0 

-4.66 

371.0 



1 (X10 






01 XJ I 







Table VI l . 

Computed 
for Case 2 

, and Quadratic Fit 

and v - 9.0. 
o 

Values 

V 

COMPUTED 

QUADRATIC 

COMPUTED 

ft/rad 

CURVE 

v = 9.0 

AIRSPEED 


(xlO A ) 

° k 

(xlO / 

ft/sec 

0.0 

0.0 

-0.13 

0.0 

1.0 

-0.25 

-0.34 

33.5 

2.0 

-0.51 

-0.57 

67.2 

3.0 

-0.82 

-0.85 

101.0 

4,0 

"1.15 

-1.17 

136,0 

5.0 

-1.52 

-1.53 

171.0 

6.0 

-1.92 

-1.92 

206.0 

7.0 

-2.36 

- 2.36 

242.0 

8.0 

-2.83 

-2.83 

279.0 

9.0 

-3.34 

-3.34 

315.0 

10.0 

-3.89 

-3-89 

351.0 

n.o 

-4.47 

-4.48 

387.0 

12.0 

*5.09 

-5.11 

422.0 

13.0 

-5.72 

-5.77 

455.0 

14.0 

-6.37 

-6.48 

487.0 

15.0 

-7.03 

-7.22 

517.0 

16.0 

-7.69 

-8.01 

546.0 



- 0,20 
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Table V 111. 

Computed ip ^ and 
for Case 1 and 

Cubic Fit Values 

v = 5.0. 
o 

• . 

V 

ft/rad 

COMPUTED 

CURVE 

(xlo\ 

QUADRATIC 
v = 5.0 

(xlO 4 ) 

COMPUTED 

AIRSPEED 

ft/sec 

0.0 

0.0 

0.0 

0 

1.0 

-0.13 

-0.13 

141 

2.0 

-0.23 

-0.21 

' 255 

3.0 

-0.29 

-0.26 

339 

4.0 

-0.28 

-0.26 

407 

5.0 

-0.21 

-0.21 

475 

6.0 

-0.13 

-0.10 

550 

7.0 

-0.06 

+0.08 

631 

8.0 

+0.01 

0.33 

724 

9.0 

0.08 

0.65 

906 

10.0 

0.15 

1.06 

925 

11.0 

0.23 

1.57 

1030 

12.0 

0.29 

2.16 

1120 

13.0 

0.35 

2.87 

1220 

14.0 

0.41 

3.68 

1310 

15.0 

0.46 

4.61 

1400 

16.0 

0.51 

5.66 

1500 



00 "0 09 '0 Oh '0 00 "0 00 "0 02*6- 

(1,-0 IX) 2 3Q0W * A 



00 

o 




t .MODE 2 (X 1 0 _tt ) 

- 0.30 0.00 0.30 0.60 0.90 1.20 
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Table IX. 

Computed ^ anc * 
for Case 2 and 

Cubic Fit Values 

v = 9.0. 

0 


V 

COMPUTED 

CUBIC FIT 

COMPUTED 

ft/ rad 

CURVE 

v - 9.0 

AIRSPEED 


(xIO 5 ) 

® c 

(xlO 5 ) 

ft/sec 

0.0 

0.0 

0.0 

0 

1.0 

-0.24 

-0.25 

161 

2,0 

-0.50 

-0.54 

3H 

3.° 

-0.77 

-0.84 

444 

4.0 

-1,03 

-1.11 

554 ' 

5.0 

-1.25 

-1.33 

641 

6.0 

-1.39 

-1.44 

708 

7.0 

-1.38 

-1.41 

759 

8.0 

-1.20 

-1.20 

800 

9.0 

-0.78 

-0.78 

832 

10.0 

-0.13 

-0.12 

864 

11.0 

+0.76 

+0.83 

898 

12.0 

1.83 

2.11 

931 

13.0 

3-03 

3.74 

966 

14.0 

4.28 

5.78 

1010 

15.0 

5.54 

8.24 

1050 

16.0 

6.77 

1.12 

1080 


•' rr’- 
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Figure 21. Quadratic Approximation of ^ for Case 1. 
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Table X. Computed and Quadratic Fit Values 
for Case 1 and v - 5*0. 


V 

ft/rad 

COMPUTED 

CURVE 

(xlO* 1 ) 

QUADRATIC 
v » 5.0 

W) 


COMPUTED 

AIRSPEED 

ft/sec 

0.0 

0.0 

+ 0.08 


0 

1.0 

-0.05 

-0.002 


161 

2.0 

-o.n 

-0.09 


327 

3.0 

-o.i8 

-0.18 


491 

4.0 

-0.27 

-0.27 


634 

5.0 

-0.36 

-0.36 


741 

6.0 

-0.45 

-0.45 


811 

7.0 

-0.54 

-0.54 


849 

8.0 

-0.65 

-0.63 


863 

9.0 

-0.92 

-0.72 


825 

10.0 

-0.96 

-0.82 


870 

n.o 

-1.14 

-0,91 


878 

12.0 

-1.34 

-1.00 


882 

13.0 

-1.54 

-1.10 


887 

14.0 

-1.75 

-1.19 


891 

15.0 

-1.98 

-1.29 


893 

16.0 

-2.21 

-1.39 


896 




CO 
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Table XI. Computed ip _ and Quadratic Fit Values , 
for Case 2"and v - 9*0. 


V 

COMPUTED 

QUADRATIC 

COMPUTED 

ft/rad 

CURVE 

v = 9.0 

AIRSPEED 


(xIO 5 ) 

o c 
(xlO b ) 

ft/sec 

0.0 

0.0 

1.04 

. 0 

1.0 

-0.08 

0.78 

180 

2.0 

-0.18 

0.50 

363 

3.0 

-0.28 

0.23 

547 

4.0 

-0.40 

-0.05 

733 

5.0 

-0.55 

-0.34 

922 

6.0 

-0.73 

-0.63 

mo 

7.0 

-0.95 

-0.91 

1300 

8.0 

-1.21 

-1.21 

1480 

9.0 

-1.50 

-1,50 

1640 

10.0 

-1.80 

-1.81 

1780 

li.o 

-2.10 

-2.11 

1890 

12.0 

-2.39 

-2.42 

1990 

13.0 

-2.66 

-2.73 

2080 

14.0 

-2.93 

-3.04 

2140 

15.0 

-3.19 

-3.36 

2200 

16.0 

-3.44 

-3.69 

2250 
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Table XU. Computed and Cubic Fit Values 

for Case 1 and v “ 5.0. 

o 


V 

ft/rad 

COMPUTED 

CURVE 

(xIO 5 ) 

CUBIC FIT 
v = 5.0 

G r 

(xlO 5 ) 

COMPUTED 

AIRSPEED 

ft/sec 

0.0 

0.0 

0.0 

0 

1.0 

-0.14 

-0.14 

362 

2.0 

-0.34 

-0.33 

690 

3.0 

-0.55 

-0.56 

964 , 

4.0 

-0.83 

-0.84 

1170 

5.0 

-1.15 

-1.15 

1300 

6.0 

-1.50 

-1.50 

1390 

7.0 

-1.87 

-1.88 

1440 

8.0 

- 2.28 

-2.29 

1480 

9.0 

-2.70 

-2.72 

1500 

10.0 

-3-16 

-3.18 

1520 

n.o 

-3.63 

-3.66 

1530 

12.0 

-4.11 

-4.15 

1530 

13.0 

-4.61 

-4.65 

1530 

14.0 

-5.12 

-5.17 

1530 

15.0 

-5.63 

-5.68 

1530 

16.0 

-6.16 

-6.21 

1530 



Ofi'O 

(s-0 
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Table XIII. Computed and Cubic Fit Values 
for Case 2 and = 9.0. 


V 

ft/ rad 

COMPUTED 

CURVE 

(xlO 6 ) 

CUBIC FIT 
= 9 .0 

0 A 

(xlO 6 ) 

COMPUTED 

AIRSPEED 

ft/sec 

0.0 

0.0 

0.0 

0 

1.0 

-0.30 

-0.18 

450 

2.0 

-0.59 

-0.43 

872 

3.0 

-0.99 

-0.68 

1230 

4.0 

-1.34 

- 0.90 

1559 

5.0 

-1.52 

-1.02 

1800 

6.0 

-1.50 

-1.01 

1990 

7.0 ' 

-1.16 

- 0.80 

2150 

8,0 . 

-0.48 

-0.34 

2310 

9.0 

+0.41 

+0.41 

2480 

10.0 

1.31 

1.51 

2660 

11.0 

2.12 

3.02 

2850 

12.0 

2.83 

4.97 

3040 

13.0 

3.43 

7.44 

3240 

14.0 

3.95 

10.45 

3440 

15.0 

4.39 

14.08 

3620 

16.0 

4.76 

18.37 

3800 
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Figure 28. t/v versus v for Case 1. 
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Figure 29. Cubic Approximation of ip r for Case 1. 
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Table XtV. 

Computed 
for Case 

ijv and Cubic Fit Values 

Kand v - 5*0. 

o 


V 

ft/rad 

COMPUTED 

CURVE 

(xIO 6 ) 

CUBIC FIT 

v = 5.0 
° A 

(xIO 6 ) 

COMPUTED 

AIRSPEED 

ft/sec 

0.0 

0.0 

0.0 

0 

1.0 

-1.37 

-1.31. 

4oo 

2.0 

-2.29 

-2.01 

755 

3.0 

-2.81 

- 2.19 

1030 

4.0 

-2.15 

-1.93 

1280 

5.0 

-1-31 

-1.31 

1550 

6.0 

-0.63 

-0.41 

1830 

7.0 

-0. 12 

+0.70 

2130 

8.0 

+0.25 

1.92 

2420 

9.0 

0.53 

3.18 

2730 

10.0 

0.72 

4.40 

3020 

n.o 

0.86 

5.50 

3320 

12.0 

0.96 

6.40 

3620 

13.0 

1.03 

7.02 

3900 

14.0 

1.08 

7.28 

4200 

15.0 

1.11 

7.09 

4500 

16.0 

1.13 

6.38 

4780 
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Table XV. 

Computed 
for Case 

and Cubic Fit Values 

2^and v = 9.0. 
o 


V 

COMPUTED 

CUBIC FIT 

COMPUTED 

ft/rad 

CURVE 

v = 9.0 

AIRSPEED 


(xlO 6 ) 

0 A 

(xlO 6 ) 

ft/sec 

0.0 

0.0 

0.0 

v 0 

1.0 

-0.18 

- 0.31 

414 

2.0 

-0.41 

-0.66 

830 

3.0 

-0.60 

-1.06 

1260 

4.0 

“0.93 

-1.52 

1680 

5.0 

-1.55 

-2.02 

2090 

6.0 

- 2.32 

-2.58 

2430 

7.0 

-3.0& 

-3.18 

2700 

8.0 

-3.81 

-3.83 

2910 

9.0 

-4.54 

-4.54 

3070 

10.0 

-5.27 

-5-28 

3200 

11.0 

-6.01 

-6.08 

3290 

12.0 

-6.78 

-6.93 

3350 

13.0 

-7.59 

-7.82 

3400 

14.0 

-8.43 

-8.75 

3430 

15.0 

-9.32 

-9.74 

3470 

16.0 

-10.30 

-10.76 

3490 
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Figure 32. ip ^ versus v for Case 1. 
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Table XVI. 

Computed 
for Case 

jjv and Cubic Fit 

l b and v = 5.0. 
o 

Val ues 

V 

ft/ rad 

COMPUTED 

CURVE 

(xlO 6 ) 

CUBIC FIT 
v - 5.0 

O r 

(xlO 6 ) 

COMPUTED 

AIRSPEED 

ft/sec 

0.0 

0.0 

0.0 

0 

1.0 

-0.63 

-0.58 

581 

2.0 

-0.88 

-0.81 

1090 

3.0 

-0.79 

-0.83 

1600 

J4.0 

-0.70 

-0.73 

2 1 60 

5.0 

-0.64 

-0.64 

2750 

6.0 

-0.60 

-0.68 

3360 

7.0 

-0.57 

-0.95 

4000 

8.0 

-0.54 

-1.59 

4650 

9.0 

-0.51 

-2.69 

5310 

10.0 

-0.49 

-4.38 

5990 

11.0 

-0.46 

-6.78 

6670 

12.0 

-0.44 

-10.00 

7360 

13.0 

-0.42 

-14.20 

8050 

14.0 

-0.40 

-19.40 

8730 

15.0 

-0.38 

-25.70 

9430 

16.0 

-0.36 

-33.40 

10100 
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Table XVII. Computed ijv and Cubic Fit Values 

for Case 2 and v = 3.0. 

o 


V 

ft/rad 

COMPUTED 

CURVE 

(xlO 7 ) 

CUBIC FIT 

v - 9.0 

° 7 

(xlO 7 ) 

COMPUTED 

AIRSPEED 

ft/sec 

0.0 

0.0 

0.0 

0 

1.0 

- 1.43 

- K39 

651 

2.0 

- 2.84 

- 2.01 

1260 

3.0 

- 4.05 

- 1.99 

1800 

4.0 

- 4.48 

- 1.49 

2250 

5.0 

- 3.14 

- 0.64 

2640 

6.0 

- 0.70 

+ 0.42 

3050 

7.0 

+ 1.25 

1.54 

3470 

8.0 

2.56 

2.59 

3910 

9.0 

3.42 

3.42 

4350 

10.0 

3.95 

3.90 

4770 

11.0 

4.26 

3.88 

5200 

12.0 

4.40 

3.22 

5630 

13-0 

4.44 

1.78 

6060 

14.0 

4.41 

- 0.58 

6490 

15.0 

4.33 

- 4.00 

6930 . 

16.0 

4.22 

- 8.62 

7360 
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ABSTRACT 


The established method for solution of the critical flutter 
velocity of the equation of motion for an aircraft structure requires 
a plot of ip versus v where ip is the imaginary part of the eigenvalue 
of the flutter equation and v is the velocity divided by the 
circular frequency of oscillation or the dependent variable of the 
flutter equation. From this plot the crossover points, where ^ is 
zero, are sought from which the lowest or critical velocity, for 
which the structure will have divergent oscillations, may be computed. 
A curve fitting approach (which is rapid, simple, and direct in 
comparison to established methods) has been developed to solve the 
flutter equation for the critical flutter velocity.’ 

The versus v curves are approximated by cubic and quadratic 
equations. The curve fitting technique utilized the first and 
second derivatives of ip with respect to v which are derived in the 
text. 

The method was tested for two structures, one structure being 
six times the total mass of the other structure. The algorithm 
never showed any tendency to diverge from the solution. The average 
time for the computation of a flutter velocity was 3-91 seconds on an 
IBM Model 50 computer for an accuracy of five per cent. For values of 
v close to the critical root of the flutter equation the algorithm 
converged on the first attempt. The maximum number, of iterations for 
convergence to the critical flutter velocity was five with an assumed 
value of v relatively distant from the actual crossover. 



